#用McCabe方法进行膜级联设计
#根据给定的VRR值进行设计
import numpy as np 
import matplotlib.pyplot as plt

Cin = np.array([10,2])#进料的A/B（染料/盐）组分浓度
Rej = np.array([0.99,0.12])
RetA = Rej[0]
RetB = Rej[1]

xf = Cin[1]/np.sum(Cin)#进料盐比例
xtarget = 0.96
xin = 1-xtarget

def fVRR(VRR):#VRR与塔板数关系
	Abta = 1-1/(1-1/VRR)*(1-(1/VRR)**RetA)
	Abtb = 1-1/(1-1/VRR)*(1-(1/VRR)**RetB)
	alpha = (1-Abtb)/(1-Abta)*(1-(1-1/VRR)*(1-Abta))/(1-(1-1/VRR)*(1-Abtb))
	xps = []#渗透侧集合
	xrs = []#渗余侧集合
	Gamas =[]
	count = 1#总膜数计数器
	countper = 0#渗透段计数器
	countret = 0#渗余段计数器
	xi = xin
	yi = xi 
	Cper = Cin
	C = Cper#迭代初值
	Ins = []
	while True:
		if xtarget <= xi:
			break
		xt = xi
		Theta = (1-1/VRR)*xi*((1-Abtb)+(1/xi-1)*(1-Abta))
		Gama = (1-Theta)/Theta 
		T = yi/(1-yi)/alpha 
		xi = T/(T+1)#平衡线的x
		plt.plot([xt,xi],[yi,yi],'b')#横线
		xrs = np.append(xrs,xi)
		Gamas = np.append(Gamas,Gama)
		In = 0
		for i in range(1,count+1):
			G = Gamas[0:i]
			In += np.sum(np.cumprod(G))
		Ins = np.append(Ins,In)
		#渗透侧操作线方程组
		def per(x,I):
			y = I/(I+1)*x+1/(1+I)*xin 
			return y 
		#渗余侧操作线方程组
		def ret(x,I):
			y = (1+1/I)*x-1/I*xtarget
			return y
		xps = np.append(xps,yi)
		I = In
		xinter = I*(1+I)/(2*I+1)*(xin/(1+I)+xtarget/I)#本级操作线交点
		#判定属于哪一段操作
		if xinter > xi:
			yi = per(xi,I)
			countper += 1
		else:
			I = Ins[count-1]
			countret += 1
			yi = ret(xi,I)
		count += 1
	y = [count,countper,countret,xrs[-1]]
	return y

VRRs = np.arange(1.1,10.01,0.01)#灵敏度分析中VRR的范围
Totals = []#总板数
xrs = []#产物在出料的质量分数
for i in VRRs:
	counts = fVRR(i)
	Totals =  np.append(Totals,counts[0])
	xrs = np.append(xrs,counts[-1])
	
xrmax = max(xrs)#找出xr最大值
position = np.where(xrs == xrmax)#寻找这个值在numpy矩阵的位置
VRRm = VRRs[position]#由于VRR与xrs一一对应，因此我们可以根据位置输出xr最大值对应的VRR
print(VRRm)
plt.subplot(1,2,1)
plt.plot(VRRs,Totals)
plt.xlabel('VRR')
plt.ylabel('Total Stages')
plt.subplot(1,2,2)
plt.plot(VRRs,xrs)
plt.xlabel('VRR')
plt.ylabel('xR')
plt.show()






	

	
	